﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Statistics;
namespace GeostatisticalTool.Classes
{
    public class TheoreticalModel:ICloneable
    {
        private Parameter msill;
        private Parameter mrange;
        private Parameter mnugget;
        private Parameter distance;
        private Func<double> regressionFunction;
        private Parameter[] regressionParameters=null;
        private Parameter[] observedParameters = null; 
        private IKrigingModel krigingModel;
        private ISpatialContinuityData spatial;
        private double[,] z;
        private bool sill=false;
        private bool range=true;
        private bool nugget=false;

        private double perc = 0.3;

        public double Perc
        {
            get { return perc; }
            set { perc = value;

            }
        }

        public event EventHandler ChangedKriginParameter;

        public ISpatialContinuityData Spatial
        {
            get { return spatial; }
            set { spatial = value;
            UpdateData();
            }
        }

        double maxd;
        double maxv;
        public void UpdateData()
        {
            double max;
            maxd   = spatial.GetDistancesValues().Max();
            maxv= spatial.GetSemivariancesValues().Max();
            max=maxd*perc;
            int ni = spatial.GetDistancesValues().Where(item => item <= max).Count();

            z = new double[2, ni];
            int nv = 0;
            for (int i = 0; i < spatial.GetDistancesValues().Length; i++)
            {
                double v = spatial.GetDistancesValues()[i];
                if (v <= max)
                {
                    z[0, nv] = v;
                    z[1, nv] = spatial.GetSemivariancesValues()[i];
                    nv++;
                }
            }
        }

        public void FitValues(bool sill, bool range, bool nugget)
        {
            this.sill=sill;
            this.range=range;
            this.nugget=nugget;
        }

        public IKrigingModel KrigingModel
        {
            get { return krigingModel; }
            set {
                  krigingModel = value;
                 OnChangedKriginParameter(EventArgs.Empty);
            }
        }

        protected virtual void OnChangedKriginParameter(EventArgs e)
        {
            if (ChangedKriginParameter != null)
                ChangedKriginParameter(this, e);
        }

        public TheoreticalModel(ISpatialContinuityData spatial)
        {
            this.spatial = spatial;
            UpdateData();
        }

        public IKrigingModel GetFunctionDefault()
        {
            GetFunction(Model.Circular, spatial.GetSemivarianceStatistic().Mean, spatial.GetLag().lagSize*spatial.GetLag().numLags/2, 0);
            return GetFunction(this.krigingModel.Range);
        }

        public IKrigingModel GetFunctionDefault(Model model)
        {
            return GetFunction(model, spatial.GetSemivarianceStatistic().Mean, spatial.GetLag().lagSize * spatial.GetLag().numLags / 2, 0);

        }

        public IKrigingModel SetFunction(Model model, double sill, double range, double nugget)
        {
            krigingModel = SelectModel(model);
            krigingModel.C0=nugget;
            krigingModel.C1=sill;
            krigingModel.Range=range;
            OnChangedKriginParameter(EventArgs.Empty);
            return krigingModel;
        }

        public static IKrigingModel SelectModel(Model value)
        {
            IKrigingModel model;
            switch (value)
            {
                case Model.Circular:
                    model = new ModelCircular();
                    break;
                case Model.Exponential:
                    model = new ModelExponential();
                    break;
                case Model.Gaussian:
                    model = new ModelGaussian();
                    break;
                case Model.Linear:
                    model = new ModelLineal();
                    break;
                case Model.Spherical:
                    model = new ModelSpherical();
                    break;
                default:
                    model = new ModelExponential();
                    break;
            }
            return model;


        }

        private void DefineParametersCalculation(){
        
            if (sill && range && nugget)
            {
                regressionParameters = new Parameter[] { mrange, msill, mnugget };
            }else if (sill && range)
            {
                regressionParameters = new Parameter[] { mrange, msill };
            } else if(sill && nugget)
            {
               regressionParameters = new Parameter[] { msill, mnugget };
            } else if (sill)
            {
             regressionParameters = new Parameter[] { msill };
            }else if (range && nugget)
            {
              regressionParameters = new Parameter[] { mrange, mnugget };
            }else if (range)
            {
              regressionParameters = new Parameter[] { mrange };
            } else if (nugget)
            {
            regressionParameters = new Parameter[] { mnugget };
            }else {
                 regressionParameters = new Parameter[] {};
           }
        }

        private void UpdateMax()
        {
            if (krigingModel.Range > maxd)
                krigingModel.Range = maxd;
            if (KrigingModel.C1 > maxv)
                KrigingModel.C1 = maxv;
        }
        private void DefineParametersCalculationValues(double sillv, double rangev, double nuggetv)
        {

        

            if (nuggetv < 0)
            {
                nuggetv = 0;
            }


            if (sill && range && nugget)
            {
               // regressionParameters = new Parameter[] { mrange, msill, mnugget };
                krigingModel.Range = regressionParameters[0].Value;
                krigingModel.C1 = regressionParameters[1].Value;
                if (regressionParameters[2].Value < 0)
                {
                    krigingModel.C0 = 0;
                }
                else
                {
                    krigingModel.C0 = regressionParameters[2].Value;
                }
                UpdateMax();
            }
            else if (sill && range)
            {
             //   regressionParameters = new Parameter[] { mrange, msill };
                krigingModel.Range = regressionParameters[0].Value;
                if (krigingModel.Range < 0)
                    krigingModel.Range = 0.1;

                krigingModel.C1 = regressionParameters[1].Value;
                krigingModel.C0 =nuggetv;
                UpdateMax();
            }
            else if (sill && nugget)
            {
              //  regressionParameters = new Parameter[] { msill, mnugget };
                krigingModel.Range = rangev;
                krigingModel.C1 = regressionParameters[0].Value;
                if (regressionParameters[1].Value < 0)
                {
                    krigingModel.C0 = 0;
                }
                else
                {
                    krigingModel.C0 = regressionParameters[1].Value;
                }
                UpdateMax();
            }
            else if (sill)
            {
              //  regressionParameters = new Parameter[] { msill };
                krigingModel.Range = rangev;
                krigingModel.C1 = regressionParameters[0].Value;
                krigingModel.C0 = nuggetv;
                UpdateMax();
            }
            else if (range && nugget)
            {
             //   regressionParameters = new Parameter[] { mrange, mnugget };
                krigingModel.Range = regressionParameters[0].Value;
                if (krigingModel.Range < 0)
                    krigingModel.Range = 0.1;

                krigingModel.C1 = sillv;
                if (regressionParameters[1].Value < 0)
                {
                    krigingModel.C0 = 0;
                }
                else
                {
                    krigingModel.C0 = regressionParameters[1].Value;
                }
                UpdateMax();

            }
            else if (range)
            {
             //   regressionParameters = new Parameter[] { mrange };
                krigingModel.Range = regressionParameters[0].Value;
                if (krigingModel.Range < 0)
                    krigingModel.Range = 0.1;
                krigingModel.C1 = sillv;
                krigingModel.C0 = nuggetv;
                UpdateMax();

            }
            else if (nugget)
            {
               // regressionParameters = new Parameter[] { mnugget };
                krigingModel.Range = rangev;
                krigingModel.C1 = sillv;
                if (regressionParameters[0].Value < 0)
                {
                    krigingModel.C0 = 0;
                }
                else
                {
                    krigingModel.C0 = regressionParameters[0].Value;
                }
                UpdateMax();
            }
            else
            {
             //   regressionParameters = new Parameter[] { };
                krigingModel.Range = rangev;
                krigingModel.C1 = sillv;
                krigingModel.C0 = nuggetv;
                UpdateMax();

            }
        }

        public IKrigingModel GetFunction(double range, double nugget)
        {
            FitValues(false,true,true);
            return GetFunction(krigingModel.ModelType, krigingModel.C1, range,nugget);
        }

        public IKrigingModel GetFunction(double range)
        {
            FitValues(true, false, false);
            return GetFunction(krigingModel.ModelType,  spatial.GetSemivarianceStatistic().Mean, range, 0);
        }


        public IKrigingModel GetFunction(Model model, double sill, double range, double nugget)
        {


             msill = new Parameter(sill);
             mrange = new Parameter(range);
             mnugget = new Parameter(nugget);
             distance = new Parameter();


            // Parameter[] regressionParameters=null;
           //  Parameter[] observedParameters = null; ;
             observedParameters = new Parameter[] { distance };
            // regressionParameters = new Parameter[] { mrange, msill,mnugget };
            DefineParametersCalculation();
            switch (model)
            { 
                case Model.Circular:
                    krigingModel = new ModelCircular();
                    regressionFunction =  () =>
                         distance > 0 && distance <= Math.Abs(this.mrange) ? this.mnugget + (((2 * Math.Abs(this.msill)) / Math.PI) * (((distance / Math.Abs(this.mrange)) * Math.Sqrt(1 - ((distance * distance) / (Math.Abs(this.mrange) * Math.Abs(this.mrange))))) + Math.Asin(distance / Math.Abs(this.mrange)))) : distance >Math.Abs(this.mrange)?this.mnugget +Math.Abs(this.msill):0
                         ;
                    break;
                case Model.Exponential:
                    krigingModel = new ModelExponential();
                    regressionFunction = () =>distance>0?
                   Math.Abs(this.mnugget) + (Math.Abs(this.msill) * (1 - Math.Exp(-3 * Math.Abs(distance) / Math.Abs(this.mrange)))):0;
                    break;
                case Model.Gaussian:
                    krigingModel = new ModelGaussian();
                    regressionFunction = () =>distance>0?
                    Math.Abs(this.mnugget) + (Math.Abs(this.msill) * (1 - Math.Exp(-3 * (Math.Abs(distance) / Math.Abs(this.mrange)) * (Math.Abs(distance) / Math.Abs(this.mrange))))):0;
                    break;
                case Model.Spherical:
                    krigingModel = new ModelSpherical();
                    regressionFunction = () =>

                                    distance > 0 && distance <= Math.Abs(this.mrange) ?
                                             Math.Abs(this.mnugget) + (Math.Abs(this.msill) * ((1.5 * (this.distance / Math.Abs(this.mrange))) -
                                             (0.5 * (this.distance / Math.Abs(this.mrange)) * (this.distance / Math.Abs(this.mrange)) * (this.distance / Math.Abs(this.mrange))))) :
                                             distance > Math.Abs(this.mrange) ?
                                                  Math.Abs(this.mnugget) + Math.Abs(this.msill) :
                                                  0; 

                    break;
            }

            if (regressionParameters.Length == 0)
            {
                krigingModel.Range =range;
                krigingModel.C1 = sill ;
                if (nugget < 0)
                {
                    krigingModel.C0 = 0;
                }
                else
                {
                    krigingModel.C0 = nugget;
                }
                OnChangedKriginParameter(EventArgs.Empty);
                return krigingModel;

            }

             LevenbergMarquardt levenbergMarquardt;

             if (regressionParameters != null)
             {
                 try
                 {

                     levenbergMarquardt = new LevenbergMarquardt(regressionFunction, regressionParameters, observedParameters, z,2);

                     for (int i = 0; i < 50; i++)
                     {
                         levenbergMarquardt.Iterate();
                     }
                     DefineParametersCalculationValues(sill, range, nugget);
                
                 }
                 catch
                 {

                     DefineParametersCalculationValues(sill, range, nugget);
                     krigingModel.Range = range;
                     krigingModel.C1 = sill;
                     if (nugget < 0)
                     {
                         krigingModel.C0 = 0;
                     }
                     else
                     {
                         krigingModel.C0 = nugget;
                     }
                     OnChangedKriginParameter(EventArgs.Empty);
                     return krigingModel;
                 
                 }

             }

             OnChangedKriginParameter(EventArgs.Empty);
             return krigingModel;

        }

        public object Clone()
        {
            TheoreticalModel model= new TheoreticalModel(this.spatial);
            model.KrigingModel = krigingModel;
            return model;
        }
    }









    public enum Model
    {
        Circular = 0,
        Exponential,
        Gaussian,
        Spherical,
        Linear
    }
}
